Surface emitter

In [1]:

# Useful for debugging
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'

import os
import yaml
import numpy as np

# Useful for debugging %load_ext autoreload %autoreload 2 %config InlineBackend.figure_format = 'retina' import os import yaml import numpy as np

In [2]:

from distgen import Generator

from distgen import Generator

In [3]:

input_yaml = """
n_particle: 100000
species: electron
start: 
  type: free
random:
  type: hammersley
total_charge:
  units: pC
  value: 10
r_dist:
  max_r:
    units: mm
    value: 1
  type: radial_uniform
p_dist:
  sigma_p:
    value: 1
    units: eV/c
  avg_p: 
    value: 100
    units: eV/c
  type: gaussian
"""

input_yaml = """ n_particle: 100000 species: electron start: type: free random: type: hammersley total_charge: units: pC value: 10 r_dist: max_r: units: mm value: 1 type: radial_uniform p_dist: sigma_p: value: 1 units: eV/c avg_p: value: 100 units: eV/c type: gaussian """

In [4]:

gen = Generator(input_yaml, verbose=True)

gen = Generator(input_yaml, verbose=True)

In [5]:

P = gen.run()

P = gen.run()

Distribution format: None
   Warning: no output file specified, defaulting to "None".
Output file: None

Creating beam distribution....
   Beam starting from: free
   Total charge: 10 pC.
   Number of macroparticles: 100000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 1 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   p distribution: Gaussian
      avg_p = 100 eV/c, sigma_p = 1.000 eV/c
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   Shifting avg_x = -4.85224E-08 mm -> 0 mm
   Scaling sigma_x = 0.499991 mm -> 0.5 mm
   Shifting avg_y = 1.44675E-06 mm -> 0 mm
   Scaling sigma_y = 0.49999 mm -> 0.5 mm
   Shifting avg_px = -3.3074E-05 eV/c -> -6.12323E-15 eV/c
   Scaling sigma_px = 57.7382 eV/c -> 57.7379 eV/c
   Shifting avg_py = 1.54491E-06 eV/c -> 0 eV/c
   Scaling sigma_py = 57.7382 eV/c -> 57.7379 eV/c
   Shifting avg_pz = -0.000118434 eV/c -> 6.12323E-15 eV/c
   Scaling sigma_pz = 57.7372 eV/c -> 57.7379 eV/c
...done. Time Elapsed: 107.992 ms.

   Created particles in .particles: 
   ParticleGroup with 100000 particles with total charge 1.0000000000000003e-11 C

In [6]:

P.plot("p")

P.plot("p")

In [7]:

P.plot("x", "px")

P.plot("x", "px")

In [8]:

P.plot("pz")

P.plot("pz")

In [9]:

P.plot("px")

P.plot("px")

In [10]:

P.plot("py")

P.plot("py")

In [11]:

from matplotlib import pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(projection="3d")

px, py, pz = P["px"], P["py"], P["pz"]

p = np.sqrt(px**2 + py**2 + pz**2)

theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

x = px / p
y = py / p
z = pz / p

ax.scatter(x[::200], y[::200], z[::200], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")

from matplotlib import pyplot as plt fig = plt.figure() ax = fig.add_subplot(projection="3d") px, py, pz = P["px"], P["py"], P["pz"] p = np.sqrt(px**2 + py**2 + pz**2) theta = np.arctan2(py, px) phi = np.arccos(pz / p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 x = px / p y = py / p z = pz / p ax.scatter(x[::200], y[::200], z[::200], ".") ax.set_xlabel(r"$\hat{p}_x$") ax.set_ylabel(r"$\hat{p}_y$") ax.set_zlabel(r"$\hat{p}_z$")

Out[11]:

Text(0.5, 0, '$\\hat{p}_z$')

In [12]:

theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

plt.plot(pcs, hist);

theta = np.arctan2(py, px) phi = np.arccos(pz / p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 plt.plot(pcs, hist);

In [13]:

rns = np.random.random(100000)

rns = np.random.random(100000)

In [14]:

pa = 0
pb = np.pi

Ca = np.cos(pa)
Cb = np.cos(pb)

pa = 0 pb = np.pi Ca = np.cos(pa) Cb = np.cos(pb)

In [15]:

rho = 1 / (Ca - Cb)

rho = 1 / (Ca - Cb)

In [16]:

phis = np.linspace(pa, pb, 1000)

phis = np.linspace(pa, pb, 1000)

In [17]:

plt.plot(phis, rho * np.sin(phis));

plt.plot(phis, rho * np.sin(phis));

In [18]:

np.trapezoid(rho * np.sin(phis), phis)

np.trapezoid(rho * np.sin(phis), phis)

Out[18]:

np.float64(0.999999175885426)

In [19]:

cdf = (Ca - np.cos(phis)) * rho

cdf = (Ca - np.cos(phis)) * rho

In [20]:

from scipy.integrate import cumulative_trapezoid as cumtrapz

from scipy.integrate import cumulative_trapezoid as cumtrapz

In [21]:

plt.plot(phis, cdf, phis, cumtrapz(rho * np.sin(phis), phis, initial=0));

plt.plot(phis, cdf, phis, cumtrapz(rho * np.sin(phis), phis, initial=0));

In [22]:

ps = np.arccos(Ca - rns * (Ca - Cb))

ps = np.arccos(Ca - rns * (Ca - Cb))

In [23]:

hist, pedges = np.histogram(ps, bins=30, density=True)

hist, pedges = np.histogram(ps, bins=30, density=True)

In [24]:

pcs = (pedges[1:] + pedges[:-1]) / 2

pcs = (pedges[1:] + pedges[:-1]) / 2

In [25]:

plt.plot(pcs, hist, phis, rho * np.sin(phis));

plt.plot(pcs, hist, phis, rho * np.sin(phis));

In [26]:

gen = Generator("data/beer.can.in.yaml", verbose=1)

# print(gen)

gen = Generator("data/beer.can.in.yaml", verbose=1) # print(gen)

In [27]:

P1 = gen.run()

P1 = gen.run()

Distribution format: gpt
Output file: beer.can.out.txt

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 10 pC.
   Number of macroparticles: 200000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 2 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   t distribution: uniform
      min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t:  1.1547 ps
   px distribution: Gaussian
      avg_px = 0 eV/c, sigma_px = 276.857 eV/c
   py distribution: Gaussian
      avg_py = 0 eV/c, sigma_py = 276.857 eV/c
   pz distribution: Gaussian
      avg_pz = 0 eV/c, sigma_pz = 276.857 eV/c
   Shifting avg_x = -2.10539E-05 mm -> 0 mm
   Scaling sigma_x = 0.999982 mm -> 1 mm
   Shifting avg_y = -1.97648E-06 mm -> 0 mm
   Scaling sigma_y = 0.999998 mm -> 1 mm
   Shifting avg_px = -0.0212516 eV/c -> 0 eV/c
   Scaling sigma_px = 276.849 eV/c -> 276.857 eV/c
   Shifting avg_py = -0.0259627 eV/c -> 0 eV/c
   Scaling sigma_py = 276.844 eV/c -> 276.857 eV/c
   Shifting avg_pz = -0.0372953 eV/c -> 0 eV/c
   Scaling sigma_pz = 276.843 eV/c -> 276.857 eV/c
   Shifting avg_t = -4.90022E-05 ps -> 0 ps
   Scaling sigma_t = 1.1547 ps -> 1.1547 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 220.904 eV/c, sigma_pz -> 166.887 eV/c
...done. Time Elapsed: 197.597 ms.

   Created particles in .particles: 
   ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C

In [28]:

P1.plot("x", "px")

P1.plot("x", "px")

In [29]:

P1.plot("p")

P1.plot("p")

In [30]:

P1.plot("pz")

P1.plot("pz")

In [31]:

P1["sigma_px"], P1["sigma_py"], P1["sigma_pz"]

P1["sigma_px"], P1["sigma_py"], P1["sigma_pz"]

Out[31]:

(np.float64(276.8570795554991),
 np.float64(276.8570795554991),
 np.float64(166.88683722575215))

In [32]:

gen = Generator("data/maxwell_boltzmann.beer.can.in.yaml", verbose=1)

gen = Generator("data/maxwell_boltzmann.beer.can.in.yaml", verbose=1)

In [33]:

P2 = gen.run()

P2 = gen.run()

Distribution format: gpt
Output file: beer.can.out.txt

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 10 pC.
   Number of macroparticles: 200000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 2 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   p distribution: Maxwell-Boltzmann
      p scale = 276.857 eV/c
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   t distribution: uniform
      min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t:  1.1547 ps
   Shifting avg_x = 1.63849E-05 mm -> 0 mm
   Scaling sigma_x = 0.999984 mm -> 1 mm
   Shifting avg_y = -2.77139E-06 mm -> 0 mm
   Scaling sigma_y = 0.999996 mm -> 1 mm
   Shifting avg_px = -0.000954618 eV/c -> -2.70524E-14 eV/c
   Scaling sigma_px = 276.853 eV/c -> 276.857 eV/rad⁰⋅⁵/c
   Shifting avg_py = 0.000856608 eV/c -> 0 eV/c
   Scaling sigma_py = 276.848 eV/c -> 276.857 eV/rad⁰⋅⁵/c
   Shifting avg_pz = -0.0142766 eV/c -> 2.70524E-14 eV/c
   Scaling sigma_pz = 276.833 eV/c -> 276.857 eV/rad⁰⋅⁵/c
   Shifting avg_t = -4.90022E-05 ps -> 0 ps
   Scaling sigma_t = 1.1547 ps -> 1.1547 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 220.909 eV/c, sigma_pz -> 166.881 eV/c
...done. Time Elapsed: 197.52 ms.

   Created particles in .particles: 
   ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C

In [34]:

P2["sigma_px"], P2["sigma_py"], P2["sigma_pz"]

P2["sigma_px"], P2["sigma_py"], P2["sigma_pz"]

Out[34]:

(np.float64(276.85707955549907),
 np.float64(276.85707955549907),
 np.float64(166.88083659251544))

In [35]:

gen = Generator("data/maxwell_boltzmann_KE.beer.can.in.yaml", verbose=1)

gen = Generator("data/maxwell_boltzmann_KE.beer.can.in.yaml", verbose=1)

In [36]:

P3 = gen.run()

P3 = gen.run()

Distribution format: gpt
Output file: beer.can.out.txt

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 10 pC.
   Number of macroparticles: 200000.
   Assuming cylindrical symmetry...

   r distribution: radial uniform
      min_r = 0 mm, max_r = 2 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   KE distribution: Maxwell-Boltzmann Energy
      kT = 150 meV
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   t distribution: uniform
      min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t:  1.1547 ps
   Shifting avg_x = 1.63849E-05 mm -> 0 mm
   Scaling sigma_x = 0.999984 mm -> 1 mm
   Shifting avg_y = -2.77139E-06 mm -> 0 mm
   Scaling sigma_y = 0.999996 mm -> 1 mm
   Shifting avg_px = -1.75487E-09 meV·s/m -> -9.02058E-20 meV·s/m
   Scaling sigma_px = 0.000922998 meV·s/m -> 0.000922972 meV·s/m
   Shifting avg_py = 1.21968E-09 meV·s/m -> 0 meV·s/m
   Scaling sigma_py = 0.000922978 meV·s/m -> 0.000922972 meV·s/m
   Shifting avg_pz = -4.48933E-08 meV·s/m -> 9.02058E-20 meV·s/m
   Scaling sigma_pz = 0.00092294 meV·s/m -> 0.000922972 meV·s/m
   Shifting avg_t = -4.90022E-05 ps -> 0 ps

   Scaling sigma_t = 1.1547 ps -> 1.1547 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 0.0007366 meV·s/m, sigma_pz -> 0.000556145 meV·s/m
...done. Time Elapsed: 205.344 ms.

   Created particles in .particles: 
   ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C

In [37]:

P3.plot("kinetic_energy")

P3.plot("kinetic_energy")

In [38]:

P3.plot("x", "px")

P3.plot("x", "px")

In [39]:

P3["sigma_px"], P3["sigma_py"], P3["sigma_pz"]

P3["sigma_px"], P3["sigma_py"], P3["sigma_pz"]

Out[39]:

(np.float64(276.70000872169317),
 np.float64(276.70000872169317),
 np.float64(166.7279863589753))

In [40]:

fig = plt.figure()
ax = fig.add_subplot(projection="3d")

px, py, pz = P3["px"], P3["py"], P3["pz"]

p = np.sqrt(px**2 + py**2 + pz**2)

theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

x = px / p
y = py / p
z = pz / p

ax.scatter(x[::500], y[::500], z[::500], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")

fig = plt.figure() ax = fig.add_subplot(projection="3d") px, py, pz = P3["px"], P3["py"], P3["pz"] p = np.sqrt(px**2 + py**2 + pz**2) theta = np.arctan2(py, px) phi = np.arccos(pz / p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 x = px / p y = py / p z = pz / p ax.scatter(x[::500], y[::500], z[::500], ".") ax.set_xlabel(r"$\hat{p}_x$") ax.set_ylabel(r"$\hat{p}_y$") ax.set_zlabel(r"$\hat{p}_z$")

Out[40]:

Text(0.5, 0, '$\\hat{p}_z$')

In [41]:

KE = np.linspace(0, 150, 10000)

E0 = 0

A1 = (0.8,)
m1 = 8

A2 = 0.1
m2 = 90

PKE = A1 * np.exp(-np.abs(KE - E0) / m1) + A2 * np.exp(-np.abs(KE - E0) / m2)
PKE = PKE / np.trapz(PKE, KE)  # Numerically intergate to normalize

plt.plot(KE, PKE)
plt.xlabel("KE (eV)")
plt.ylabel("$\\rho(KE)$ (1/eV)");

KE = np.linspace(0, 150, 10000) E0 = 0 A1 = (0.8,) m1 = 8 A2 = 0.1 m2 = 90 PKE = A1 * np.exp(-np.abs(KE - E0) / m1) + A2 * np.exp(-np.abs(KE - E0) / m2) PKE = PKE / np.trapz(PKE, KE) # Numerically intergate to normalize plt.plot(KE, PKE) plt.xlabel("KE (eV)") plt.ylabel("$\\rho(KE)$ (1/eV)");

/tmp/ipykernel_2270/3289071326.py:12: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
  PKE = PKE / np.trapz(PKE, KE)  # Numerically intergate to normalize

In [42]:

dat = np.zeros((len(KE), 2))
dat[:, 0], dat[:, 1] = KE, PKE

np.savetxt("KEdist.txt", dat, header="KE     PKE", comments="")

dat = np.zeros((len(KE), 2)) dat[:, 0], dat[:, 1] = KE, PKE np.savetxt("KEdist.txt", dat, header="KE PKE", comments="")

In [43]:

q = 1000 * 1.60217663e-19

print(q)

input_yaml = """
n_particle: 1000
species: electron
start: 
  type: cathode
random:
  type: hammersley
total_charge:
  units: C
  value: 1.60217663e-17
r_dist:
  max_r:
    units: mm
    value: 14.6
  type: radial_uniform
KE_dist:
  file: KEdist.txt
  units: eV
  type: file1d
transforms:
  sx:
    avg_x:
      units: millimeter
      value: -50
    type: set_avg x
  sy:
    avg_y:
      units: millimeter
      value: 50
    type: set_avg y
  sz:
    avg_z:
      units: millimeter
      value: 65
    type: set_avg z
output:
    file: None
"""

q = 1000 * 1.60217663e-19 print(q) input_yaml = """ n_particle: 1000 species: electron start: type: cathode random: type: hammersley total_charge: units: C value: 1.60217663e-17 r_dist: max_r: units: mm value: 14.6 type: radial_uniform KE_dist: file: KEdist.txt units: eV type: file1d transforms: sx: avg_x: units: millimeter value: -50 type: set_avg x sy: avg_y: units: millimeter value: 50 type: set_avg y sz: avg_z: units: millimeter value: 65 type: set_avg z output: file: None """

1.6021766299999998e-16

In [44]:

inputs = yaml.safe_load(input_yaml)

n = 1000

inputs["n_particle"] = n
inputs["output"]["file"] = "test_ion_writer.ion"
inputs["output"]["type"] = "simion"
inputs["total_charge"]["value"] = n * 1.60217663e-19
inputs["transforms"]["sx"]["avg_x"]["value"] = -50
inputs["transforms"]["sy"]["avg_y"]["value"] = 50
inputs["transforms"]["sz"]["avg_z"]["value"] = 55

inputs = yaml.safe_load(input_yaml) n = 1000 inputs["n_particle"] = n inputs["output"]["file"] = "test_ion_writer.ion" inputs["output"]["type"] = "simion" inputs["total_charge"]["value"] = n * 1.60217663e-19 inputs["transforms"]["sx"]["avg_x"]["value"] = -50 inputs["transforms"]["sy"]["avg_y"]["value"] = 50 inputs["transforms"]["sz"]["avg_z"]["value"] = 55

In [45]:

gen = Generator(inputs, verbose=1)

gen = Generator(inputs, verbose=1)

In [46]:

P = gen.run()

P = gen.run()

Distribution format: simion
Output file: /home/runner/work/distgen/distgen/docs/examples/test_ion_writer.ion

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 1.60218E-16 C.
   Number of macroparticles: 1000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 14.6 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   KE distribution: KE-distribution file: "/home/runner/work/distgen/distgen/docs/examples/KEdist.txt"
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   Shifting avg_x = 0.00564031 mm -> 0 mm
   Scaling sigma_x = 7.29551 mm -> 7.3 mm
   Shifting avg_y = 0.00120927 mm -> 0 mm
   Scaling sigma_y = 7.29741 mm -> 7.3 mm
   Shifting avg_px = -6.09774E-09 eV·s/m -> -9.89015E-22 eV·s/m
   Scaling sigma_px = 1.11822E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_py = -2.50205E-09 eV·s/m -> 0 eV·s/m
   Scaling sigma_py = 1.11987E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_pz = -2.52855E-08 eV·s/m -> 9.89015E-22 eV·s/m
   Scaling sigma_pz = 1.11509E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 8.07962E-06 eV·s/m, sigma_pz -> 7.72348E-06 eV·s/m
   Applying user supplied transform: "sx" = set_avg x...
      Setting avg_x -> -50 mm.
   Applying user supplied transform: "sy" = set_avg y...
      Setting avg_y -> 50 mm.
   Applying user supplied transform: "sz" = set_avg z...
      Setting avg_z -> 55 mm.
...done. Time Elapsed: 26.6883 ms.

   Created particles in .particles: 
   ParticleGroup with 1000 particles with total charge 1.6021766300000003e-16 C

In [47]:

P.plot("x", "y", figsize=(5, 5))

P.plot("x", "y", figsize=(5, 5))

In [48]:

P.plot("x", "px", figsize=(5, 5))

P.plot("x", "px", figsize=(5, 5))

In [49]:

P.plot("kinetic_energy")

P.plot("kinetic_energy")

In [50]:

from distgen.writers import writer

writer("simion", gen.beam(), "test_ion_writer.ion", params={"color": 0})

from distgen.writers import writer writer("simion", gen.beam(), "test_ion_writer.ion", params={"color": 0})

Distribution format: simion
Output file: /home/runner/work/distgen/distgen/docs/examples/test_ion_writer.ion

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 1.60218E-16 C.
   Number of macroparticles: 1000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 14.6 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   KE distribution: KE-distribution file: "/home/runner/work/distgen/distgen/docs/examples/KEdist.txt"
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   Shifting avg_x = 0.00564031 mm -> 0 mm
   Scaling sigma_x = 7.29551 mm -> 7.3 mm
   Shifting avg_y = 0.00120927 mm -> 0 mm
   Scaling sigma_y = 7.29741 mm -> 7.3 mm
   Shifting avg_px = -6.09774E-09 eV·s/m -> -9.89015E-22 eV·s/m
   Scaling sigma_px = 1.11822E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_py = -2.50205E-09 eV·s/m -> 0 eV·s/m
   Scaling sigma_py = 1.11987E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_pz = -2.52855E-08 eV·s/m -> 9.89015E-22 eV·s/m
   Scaling sigma_pz = 1.11509E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 8.07962E-06 eV·s/m, sigma_pz -> 7.72348E-06 eV·s/m
   Applying user supplied transform: "sx" = set_avg x...
      Setting avg_x -> -50 mm.
   Applying user supplied transform: "sy" = set_avg y...
      Setting avg_y -> 50 mm.
   Applying user supplied transform: "sz" = set_avg z...
      Setting avg_z -> 55 mm.
...done. Time Elapsed: 28.3866 ms.

In [51]:

from scipy.constants import physical_constants

mc2 = 1e6 * physical_constants["electron mass energy equivalent in MeV"][0]
e_ = physical_constants["elementary charge"][0]
me = physical_constants["electron mass in u"][0]


def particle_group_to_SIMION(P, filename, color=0):
    header = ";0"

    simion_params = [
        "TOB",
        "MASS",
        "CHARGE",
        "X",
        "Y",
        "Z",
        "AZ",
        "EL",
        "KE",
        "CWF",
        "COLOR",
    ]

    # simion_units = {
    #    "TOB": "usec",
    #    "MASS": "amu",
    #    "CHARGE": "e",
    #    "X": "mm",
    #    "Y": "mm",
    #    "Z": "mm",
    #    "AZ": "deg",
    #    "EL": "deg",
    #    "CWF": "",
    #    "COLOR": "",
    # }

    data = np.zeros((len(P), len(simion_params)))

    data[:, simion_params.index("TOB")] = P.t * 1e6  # [P.t] = sec, convert to usec

    if P.species == "electron":
        data[:, simion_params.index("MASS")] = np.full(len(P), me)
        data[:, simion_params.index("CHARGE")] = np.full(len(P), -1)
    else:
        raise ValueError(f"Species {P.species} is not supported")

    data[:, simion_params.index("X")] = P.z * 1e3
    data[:, simion_params.index("Y")] = P.y * 1e3
    data[:, simion_params.index("Z")] = -P.x * 1e3

    px = P.pz
    py = P.py
    pz = -P.px

    data[:, simion_params.index("KE")] = P.kinetic_energy  # [eV]
    data[:, simion_params.index("AZ")] = np.arctan2(-pz, px) * (180 / np.pi)  # [deg]
    data[:, simion_params.index("EL")] = np.arctan2(py, np.sqrt(px**2 + pz**2)) * (
        180 / np.pi
    )  # [deg]

    data[:, simion_params.index("CWF")] = (
        P.weight / e_
    )  # Charge Weighting Factor, derive from particle group weights
    data[:, simion_params.index("COLOR")] = np.full(len(P), color)
    # fname, X, fmt='%.18e', delimiter=' '

    np.savetxt(filename, data, delimiter=",", header=header, comments="", fmt="  %.9e")

from scipy.constants import physical_constants mc2 = 1e6 * physical_constants["electron mass energy equivalent in MeV"][0] e_ = physical_constants["elementary charge"][0] me = physical_constants["electron mass in u"][0] def particle_group_to_SIMION(P, filename, color=0): header = ";0" simion_params = [ "TOB", "MASS", "CHARGE", "X", "Y", "Z", "AZ", "EL", "KE", "CWF", "COLOR", ] # simion_units = { # "TOB": "usec", # "MASS": "amu", # "CHARGE": "e", # "X": "mm", # "Y": "mm", # "Z": "mm", # "AZ": "deg", # "EL": "deg", # "CWF": "", # "COLOR": "", # } data = np.zeros((len(P), len(simion_params))) data[:, simion_params.index("TOB")] = P.t * 1e6 # [P.t] = sec, convert to usec if P.species == "electron": data[:, simion_params.index("MASS")] = np.full(len(P), me) data[:, simion_params.index("CHARGE")] = np.full(len(P), -1) else: raise ValueError(f"Species {P.species} is not supported") data[:, simion_params.index("X")] = P.z * 1e3 data[:, simion_params.index("Y")] = P.y * 1e3 data[:, simion_params.index("Z")] = -P.x * 1e3 px = P.pz py = P.py pz = -P.px data[:, simion_params.index("KE")] = P.kinetic_energy # [eV] data[:, simion_params.index("AZ")] = np.arctan2(-pz, px) * (180 / np.pi) # [deg] data[:, simion_params.index("EL")] = np.arctan2(py, np.sqrt(px**2 + pz**2)) * ( 180 / np.pi ) # [deg] data[:, simion_params.index("CWF")] = ( P.weight / e_ ) # Charge Weighting Factor, derive from particle group weights data[:, simion_params.index("COLOR")] = np.full(len(P), color) # fname, X, fmt='%.18e', delimiter=' ' np.savetxt(filename, data, delimiter=",", header=header, comments="", fmt=" %.9e")

In [52]:

particle_group_to_SIMION(P, "text_ion_file.ion")

particle_group_to_SIMION(P, "text_ion_file.ion")

In [53]:

def read_simion_ION_file(filename):
    data = np.loadtxt(filename, comments=";", delimiter=",", skiprows=1)

    simion_params = [
        "TOB",
        "MASS",
        "CHARGE",
        "X",
        "Y",
        "Z",
        "AZ",
        "EL",
        "KE",
        "CWF",
        "COLOR",
    ]

    return {simion_params[ii]: data[:, ii] for ii, p in enumerate(simion_params)}

def read_simion_ION_file(filename): data = np.loadtxt(filename, comments=";", delimiter=",", skiprows=1) simion_params = [ "TOB", "MASS", "CHARGE", "X", "Y", "Z", "AZ", "EL", "KE", "CWF", "COLOR", ] return {simion_params[ii]: data[:, ii] for ii, p in enumerate(simion_params)}

In [54]:

ions1 = read_simion_ION_file("text_ion_file.ion")
ions2 = read_simion_ION_file("test_ion_writer.ion")

ions1 = read_simion_ION_file("text_ion_file.ion") ions2 = read_simion_ION_file("test_ion_writer.ion")

In [55]:

fig = plt.figure()
ax = fig.add_subplot(projection="3d")

px, py, pz = P["px"], P["py"], P["pz"]

p = np.sqrt(px**2 + py**2 + pz**2)

theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

x = px / p
y = py / p
z = pz / p

ax.scatter(x[::], y[::], z[::], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")

fig = plt.figure() ax = fig.add_subplot(projection="3d") px, py, pz = P["px"], P["py"], P["pz"] p = np.sqrt(px**2 + py**2 + pz**2) theta = np.arctan2(py, px) phi = np.arccos(pz / p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 x = px / p y = py / p z = pz / p ax.scatter(x[::], y[::], z[::], ".") ax.set_xlabel(r"$\hat{p}_x$") ax.set_ylabel(r"$\hat{p}_y$") ax.set_zlabel(r"$\hat{p}_z$")

Out[55]:

Text(0.5, 0, '$\\hat{p}_z$')

In [56]:

for p in ions1:
    print(f"{p}:", max(np.abs(ions1[p] - ions2[p])))

for p in ions1: print(f"{p}:", max(np.abs(ions1[p] - ions2[p])))

TOB: 0.0
MASS: 0.0
CHARGE: 0.0
X: 0.0
Y: 0.0
Z: 0.0
AZ: 0.0
EL: 0.0
KE: 0.0
CWF: 0.0
COLOR: 0.0

In [57]:

plt.plot(ions1["KE"], ions1["KE"] / ions2["KE"], ".")
plt.xlabel("KE (eV)");

plt.plot(ions1["KE"], ions1["KE"] / ions2["KE"], ".") plt.xlabel("KE (eV)");

In [58]:

os.remove("text_ion_file.ion")
os.remove("KEdist.txt")
os.remove("test_ion_writer.ion")

os.remove("text_ion_file.ion") os.remove("KEdist.txt") os.remove("test_ion_writer.ion")

In [ ]:

In [ ]: